A study of the influence of isotopic substitution on the melting point 
and temperature of maximum density of water by means of path inte- 
gral simulations of rigid models. 
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The melting point of ice Ih, as well as the temperature of maximum density (TMD) in the liquid phase, has been computed using 
the path integral Monte Carlo method. Two new models are introduced; TIP4PQ_D20 and TIP4PQ_T20 which are specifically 
designed to study D2O and T2O respectively. We have also used these models to study the "competing quantum effects" proposal 
of Habershon, Markland and Manolopoulos; the TIP4PQ/2005, TIP4PQ/2005 (D 2 0) and TIP4PQ/2005 (T 2 0) models are able to 
study the isotopic substitution of hydrogen for deuterium or tritium whilst constraining the geometry, while the TIP4PQ_D20 and 
TIP4PQ_T20 models, where the O-H bond lengths are progressively shortened, permit the study of the influence of geometry 
(and thus dipole moment) on the isotopic effects. For TIP4PQ_D20 - TIP4PQ/2005 we found a melting point shift of 4.9 K 
(experimentally the value is 3.68K) and a TMD shift of 6K (experimentally 7.2K). For TIP4PQ_T20 - TIP4PQ/2005 we found a 
melting point shift of 5.2 K (experimentally the value is 4.49K) and a TMD shift of 7K (experimentally 9.4K). 



1 Introduction 

Water molecules are composed of two hydrogen atoms and one 
oxygen atom (H 2 0). One can substitute the hydrogen atom for 
isotopes of hydrogen: one can replace hydrogen for deuterium 
to form deuterium oxide ( 2 H 2 0, or D 2 0), popularly known as 
"heavy water", or with tritium, forming tritium oxide ( 3 H 2 0, 
or T 2 0). These isotopically substituted forms of water differ 
in their physical properties, such as heat capacities (C p ), melt- 
ing point (T m ), diffusion coefficients and the temperature of the 
maximum in density (TMD). 

Within the Born-Oppenheimer approximation the "adiabatic 
surface", or potential energy surface (PES) is unaffected by 
such isotopic substitutions; any shift in experimental properties 
is due to a different probability distribution of the different iso- 
topes on the same PES, something that classical statistical me- 
chanics is unable to describe. The quantum nature of the nuclei 
becomes increasingly relevant when dealing with light atoms, 
in particular hydrogen. For this reason the incorporation of nu- 
clear quantum effects in water is germane. A complete quan- 
tum mechanical description of water would require the solution 
of the electronic Schroedinger equation for the electronic part 
of the Born-Oppenheimer approximation, in conjunction with, 
say, the path-integral formalism for the nuclear contribution—"—. 
Such a "complete" solution is still in the far future 4 . An inter- 
mediate approach is to use an empirical potential in place of 
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the potential energy surface. A great many empirical potentials 
exist for water 5 (perhaps more than for any other molecule), 
having varying degrees of success 6 . One recently proposed 
classical model has been shown to be capable of reproducing 
a good number of the thermodynamic and transport properties 
of water, namely the TIP4P/2005 modeli 

Since in water nuclear quantum effects are significant 8 ' 9 , one 
must conclude that the parameters of these classical models im- 
plicitly include these quantum contributions. A path integral 
simulation of a model such as TIP4P/2005 would be inappro- 
priate as it would lead to a "double-counting" of the nuclear 
quantum effects, so in view of this a variant of this model was 
developed; namely the TIP4PQ/2005 model 10 . It was found 
that an increase of 0.02<? in the charge of the proton led to one 
of the most quantitative phase diagrams of water calculated to 
dateii. 

Recently a very interesting suggestion has been put forward 
by Habershon, Markland and Manolopoulos 12 , that of "com- 
peting quantum effects" in water; they have proposed that zero 
point fluctuations lead to a longer O-H bond length, and thus a 
larger dipole moment, making the water molecule "less" quan- 
tum, whereas on the other hand inter-molecular quantum fluc- 
tuations serve to weaken the hydrogen bonds, making the liquid 
as a whole more "quantum". An analogous process almost cer- 
tainly takes place upon isotopic substitution; the replacement of 
hydrogen with deuterium has two effects: on the one hand the 
hydrogen bond becomes stronger, since D is less delocalised, 
i.e. more classical, than H, whilst on the other hand replacing 
H with D reduces the intramolecular OH covalent bond length, 
which in turn decreases the dipole moment of the molecule, 
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2. 1 Calculation of the melting point 



2 METHODOLOGY AND SIMULATION DETAILS 



effectively reducing the strength of the hydrogen bond. 

To study these competing quantum effects, and to exam- 
ine the influence of isotopic substitution on both the melting 
point of ice Ih and the location of the temperature of maximum 
density (TMD) in the liquid, we studied the aforementioned 
TIP4PQ/2005 model along with two new models, specifically 
designed for simulations of D?0 and T2O. Each of these mod- 
els are both rigid and non-polarisable. Water molecules are, 
beyond a doubt, flexible in nature. Furthermore, it is known 
that the isotopic substitutions noticeably affect the vibrational 
properties of water 13 and also that there is a degree of cou- 
pling between intermolecular and intramolecular modes. Since 
intramolecular vibrations are generally high frequency oscilla- 
tions a quantum rather a classical description of these vibration 
would be desirable 14 . That said however, the thermodynamic 
properties of the condensed phases of water are largely domi- 
nated by the intermolecular hydrogen bond rather than by the 
intramolecular vibrations. For this reason an analysis of how 
far one can go in the description of isotopic effects on the TMD 
and melting point of water using rigid models is still pertinent. 
It will be shown that when the same model is used to study 
each of the isotopes of water the variation of the properties, 
although qualitatively correct, are overestimated. However, it 
will be shown that by simply shortening the O-D and O-T bond 
length with respect to that of O-H, then predictions of isotopic 
effects are in reasonable agreement with experimental results. 

2 Methodology and simulation details 

Recent experiments have indicated that the O-D bond length is 
shorter than the O-H bond length by w 0.5% 15 . In view of this 
we have taken the TIP4PQ/2005 model, and also shortened the 
O-D distance by a similar amount (by 0.004A to be precise), 
resulting in the TIP4PQ_D20 model. The location of the neg- 
ative charge (situated on the massless site M) was also shifted 
so as to maintain the same ratio of afoivi/^OH as in the original 
TIP4PQ/2005 model. On doing this the relative distances be- 
tween the charges (responsible of the hydrogen bond strength) 
and the Lennard-Jones site (which controls the short range re- 
pulsive forces) remains unchanged, and provides a dipole to 
quadrupole ratio close to one, which has been shown to lead to 
a good phase diagram for TIP4P-like models 16 . We have also 
parameterised a model for T2O (TIP4PQ_T20) along the same 
lines. Given the paucity of experimental data for the O-T bond 
length in the liquid phase, we have taken the liberty of short- 
ening the O-T bond length by 0.006A with respect to the H-O 
bond length, again maintaining the bond length ratio c/om/^oh- 
The resulting parameters are given in TableQ] It is worth reiter- 
ating that all of these models are rigid and non-polarisable. The 
path integral methodology for rigid rotors was employed and 
has been discussed in detail elsewhere 17 , and we shall restrict 



Model 


A (cm l ) 


B (cm : ) 


C (cm l ) 


TIP4PQ/2005 


27.432 


14.595 


9.526 


TIP4PQ/2005 (D 2 0) 


15.262 


7.303 


4.939 


TIP4PQ/2005 (T 2 0) 


11.211 


4.877 


3.398 


TIP4PQ_D20 


15.390 


7.365 


4.981 


TIP4PQ_T20 


11.353 


4.939 


3.441 



Table 2 Rotational constants for each of the models. 



ourselves to describing the most salient aspects of the work un- 
dertaken here. The NVM propagator 18 , exact for asymmetric 
tops, was used. The NVM propagator is based on the work of 
Muser and Berne for symmetric tops 19 . We used P = 7 Trotter 
slices, or "replicas", for all simulations. In Table [2] the rota- 
tional constants for the various models are presented. Simula- 
tions of the liquid phase consisted of 360 molecules, and the 
ice Ij, phase consisted of 432 molecules. The proton disordered 
configuration of ice \, having both zero dipole moment as well 
as satisfying the Bernal-Fowler rules 20 was obtained by means 
of the algorithm of Buch et al.—. The Lennard-Jones part of 
the potential was truncated at 8. 5 A and long range corrections 
were added. Coulombic interactions were treated using Ewald 
summation method. In the ice Ij, phase the NpT ensemble was 
used, with anisotropic changes in the volume of the simulation 
box; each side being able to fluctuate independently. All sim- 
ulations were performed at a pressure of 1 bar. A Monte Carlo 
cycle consists of a trial move per particle (the number of parti- 
cles is equal to NP where N is the number of water molecules) 
plus a trial volume change in the case of NpT simulations. 

2.1 Calculation of the melting point 

Calculation of the melting point of ice Ih consists of three steps: 
Step 1: The first step is to calculate the classical melting 
point for the model of interest. This consists in calculating 
the free energy of the solid phase via Einstein crystal calcu- 
lations, followed by the addition of the residual entropy as cal- 
culated by Pauling 22 . To obtain the free energy of the fluid 
phase a thermodynamic path is constructed, making a connec- 
tion to a Lennard-Jones reference fluid whose free energy is 
well known. Once the free energies of the fluid and solid phases 
of the classical system are known for a reference state thermo- 
dynamic integration is performed to obtain the temperature for 
which both phases have the same chemical potential (at stan- 
dard pressure), i.e. the classical melting point. A thorough 
description of this procedure can be found in Ref. 23 . 

Step 2: At this classical melting point one then calculates the 
difference in chemical potential between ice Ih and water for 
the quantum system (A/i) via: 

IeT = L I 7 {m^f) ~ { m^r) dA (1) 
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3 RESULTS 



2.2 Calculation of the TMD 



Table 1 Parameters for TIP4PQ/2005 and the new TIP4PQJD20 and TIP4PQ_T20 models. The distance between the oxygen and hydrogen 
sites is don. The angle, in degrees, formed by hydrogen, oxygen, and the other hydrogen atom is denoted by H-O-H. The Lennard- Jones site 
is located on the oxygen with parameters o~ and e. The charge on the proton is q^. The negative charge is placed in a point M at a distance 
dowi from the oxygen along the H-O-H bisector. 



Model 


d H (A) 


ZH-O-H 


cr(A) 


e/k B (K) 


<?H(e) 


^om(A) 


p (Debye) 


TIP4PQ/2005 


0.9572 


104.52 


3.1589 


93.2 


0.5764 


0.1546 


2.388 


TIP4PQ_D20 


0.9532 


104.52 


3.1589 


93.2 


0.5764 


0.153954 


2.378 


TIP4PQ.T20 


0.9512 


104.52 


3.1589 


93.2 


0.5764 


0.1536 


2.373 



where K represents the total kinetic energy, given by: 

K = K tm + K rot , (2) 



where 



Kf, 



3NP 



MP 



N P 

LEW 



(3) 



NpT 



where P is the number of 'beads" and M is the molecular mass, 
and 



K„ 



1 p N 1 

p E H e rot,i 
r r=l;=l , 



(4) 



NpT 



where e^ 1 is the rotational energy term of the NVM propaga- 
tor (for details see Ref.—). The parameter A' is defined so that 
the mass of each atom i of the system is scaled as m; = nij $/X' 
where m,-,o is the mass of atom i in the original system. The 
values of m^o for H, D, T and O were taken from Ref.—. On 
increasing the atomic masses by the factor 1/A' the geometry 
and centre of mass of the molecule remains unchanged. Simi- 
larly the eigenvalues of the inertia tensor, and thus the energies 
of the asymmetric top appearing in the rotational propagator, 
are also scaled by this factor. Such a linear scaling particu- 
larly convenient for practical purposes. However, the same is 
not true for a transformation from, say, the TIP4PQ/2005 to 
TIP4PQ_D20 models, since the geometry and mass distribu- 
tion varies between the models. Similarly there is no simple 
scaling for the values of the rotational constants A, B and C 
used in the calculation of the propagator. For this reason we 
perform the integration to infinite mass for each of the models, 
rather than perform a transformation between models. Eq. [TJ 
embodies the idea that the phase that has the higher kinetic en- 
ergy will also have the higher chemical potential and as a result 
will become less stable in the quantum system. Note that for 
the TIP4PQ/2005 model the melting point of the classical sys- 
tem is the same for H2O, D2O and T2O since the melting point 
of a classical system is independent of the molecular mass. The 
integrand of Eq. QJrepresents the transformation from H2O (or 
D2O or T2O) to an infinitely massive molecule of water. This 
integral was evaluated using seven values of A' between 1/7 and 



1 (i.e. A' = 1,6/7,5/7,4/7,3/7,2/7, and 1/7) by performing 
runs of about one million cycles at each value of A'. We do not 
go beyond A' = 1/7 due to the increased expense in the eval- 
uation of the propagator, however, in Ref.— it is shown that 
the integral is well behaved down to A = for the case of the 
harmonic oscillator. Furthermore our direct coexistence simu- 
lation s U ' 25 corroborate our melting point calculations, which 
take to be an indication that there is no "anomalous" behaviour 
in the region between A' = 1/7 and A' = 0. 

Step 3: Again using thermodynamic integration 23 , the free 
energy of each phase of the quantum system is determined as a 
function of T: 



G(T 2 ,p) _ G(T\,p) 
Nk B T 2 Nk B Ti 



f T 2 H(T) 
J Tl Nk B T 2 



dT 



(5) 



where G is the Gibbs energy function and H is the enthalpy. 
This provides the location of the melting point of the quantum 
system as the temperature at which the chemical potential of 
ice Ih and water become identical. 

2.2 Calculation of the TMD 

The determination of the location of the TMD consisted in par- 
ticularly long simulation runs (up to 9 million Monte Carlo cy- 
cles per temperature) for a range of temperatures that bracket 
the location of the TMD. Once the densities as a function of 
temperature were obtained, they were fitted to a quadratic poly- 
nomial, whose maxima was taken to be the location of the 
TMD. 



3 Results 

3.1 Melting point of the TIP4PQ/20005 model 

The classical value of the melting point for the TIP4PQ/20005 
was calculated to be 282K 1 1 . As per Step 2 of the methodology 
outlined previously, the integrand of equationQJwas calculated 
and the results are presented in Figure [TJ 

One can see that this integrand is positive, indicating that the 
molecules have more kinetic energy in the ice phase than they 
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3.2 Isotope effects on the melting point 



3 RESULTS 



Model T(K) Aju(I h - liquid) / (k B T) 



TIP4PQ/2005 


282 


0.198 


TIP4PQ/2005 (D 2 0) 


282 


0.120 


TIP4PQ/2005 (T 2 0) 


282 


0.092 


TIP4PQ_D20 


276 


0.108 


TIP4PQ_T20 


273 


0.084 



Table 3 The difference in the chemical potential between the ice and liquid phases in the quantum system, evaluated at the T m of the classical 
system at a pressure of p = 1 bar. (Error ± 0.01). 



Model 


T 




?tmd 


i TMD I TMD 


H 2 (experiment) 


273.15 





277.13 





D2O (experiment) 


276.83 


3.68 


284.34 


7.2 


T2O (experiment) 


277.64 


4.49 


286.55 


9.4 






T TIP4PQ/20U5 


?tmd 


_ rr TIP4PQ/2()05 
I TMD - 1 TMD 


TIP4PQ/2005 


258.3 





284 





TIP4PQ/2005 (D 2 0) 


267.7 


9.4 


295 


11 


TIP4PQ/2005 (T 2 0) 


271.8 


13.5 


300 


16 


TIP4PQ/2005 


258.3 





284 





TIP4PQ_D20 


263.2 


4.9 


290 


6 


TIP4PQ_T20 


263.5 


5.2 


291 


7 



Table 4 Melting points and temperatures of maximum density of the models (all temperatures are in Kelvin). 



do in the liquid phase, thus ice I], is less stable in the quan- 
tum system, which in turn implies that the melting point will 
move to lower temperatures in the quantum system. Since the 
hydrogen bonds are stronger in the ice phase quantum effects 
are more influential in the ice phase. The integrand is fairly 
smooth, and so can be fitted to a straight line down to very small 
values of A', and it is from this fit that we obtain the value of the 
integral. The values for these integrals are presented in Table 
[3] Having these integrals we proceeded to Step 3, i.e. the ther- 
modynamic integration given in Eq. [5] To do this path integral 
simulations were performed for both ice Ih and water at vari- 
ous temperatures along the p = 1 bar isobar. The melting point 
of the quantum system is the temperature at which the chemi- 
cal potential of both I], and water are the same. The resulting 
melting points are given in Table [4] 

The melting point of the TIP4PQ/2005 model is 258K, 
which is approximately 15K below the experimental value. 
TIP4PQ/2005 is not alone in underestimating the melting 
point; the flexible q-TIP4P/F model—, also designed for use 
in path integral simulations, has a similar melting point (251 
K— ). This is probably due to the fact that both the q- 
TIP4P/F and TIP4PQ/2005 models are derived from the classi- 
cal TIP4P/2005 model which has T,„ = 252 K. The TTM2.1-F 
and TTM3-F models 27 , which are both flexible and polaris- 
able and were obtained from fits to high level ab initio calcu- 
lations, have somewhat lower melting points; 228 K 28 and 225 
K— respectively, while the q-SPC/Fw model 30 has a T m of 195 



K— . Conversely, density functional theory predictions for the 
melting point tend to significantly overestimate the experimen- 
tal value; two common functionals, PBE0 and BLY3P, 31 have 
a melting point of T„, =415 K. 

From Table |5]one can see that the TIP4PQ/2005 models un- 
derestimate the melting enthalpy (1.099 kcal/mol, whereas the 
experimental value is 1.436 kcal/mol). The enthalpy of melting 
was obtained from NpT simulations of both the solid phase and 
the liquid phase, for each model, both at the melting point, then 
simply taking the enthalpy difference at this temperature. Both 
of these results, the melting point and the melting enthalpy, 
were also underestimated in classical simulations of the classi- 
cal model TIP4P/2005 (251 K and 1.15 kcal/mol respectively). 
From this we can deduce that the inclusion of nuclear quantum 
effects has relatively little influence over these properties, and 
that any discrepancy with experiment is due to the approximate 
description of the PES implied in the empirical TIP4P/2005 and 
TIP4PQ/2005 models. 

3.2 Isotope effects on the melting point 

In classical simulations the melting point is independent of the 
molecular mass, thus the melting points of the TIP4PQ/2005 
(D 2 0) and TIP4PQ/2005 (T 2 0) models is the same as that of 
the TIP4PQ/2005 model, namely 282K. As per Step 2 of the 
methodology outlined previously the integrand of equation Q] 
was calculated (see Figure [TJ and the integral of Eq.(l) eval- 
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3 RESULTS 



3.4 The temperature difference between the melting point and the TMD 
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0.2 
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0.1 
0.05 




TIP4PQ/2005 




TIP4PQ D20 



0.2 



0.4 



0.6 



0.8 



X 



Figure 1 Integrand of Eq. Q](i.e. (A" Ih - Ari iquid ) /(X'Nk B T)) as a 
function of A'. The integral of the curves (from to 1) yields A/i(Ih • 
liquid)/ (k B T). Key: TIP4PQ/2005 red line, TIP4PQJD20 blue 
dashed line. 
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Figure 2 Plot of the isobars (p — 1 bar) used to calculate the TMD 
(points), along with fits to experimental results 37 for H2O and D2O 
(lines) (the authors were unable to locate experimental data for T2O). 



uated (see Table [3]). Thermodynamic integration was then un- 
dertaken leading to the melting points of the quantum system; 
268K for TIP4PQ/2005 (D 2 0) and 272K for TIP4PQ/2005 
(T2O). This increase in the melting point qualitatively mirrors 
experimental results, however, the magnitude of the shift is over 
estimated (see column 3 of Table 4). 

The same procedure was applied to the TIP4PQ_D20 and 
TIP4PQ_T20 models, which have classical melting points of 
276 K and 273 K respectively. We find a difference of 4.9 
K between TIP4PQJ320 and TIP4PQ/2005, and 5.2 K be- 
tween TIP4PQ.T20 and TIP4PQ/2005, which is in much bet- 
ter agreement with the experimental value of the shift. Similar 
values for the melting point differences were found for the q- 
TIP4P/F model; 6.5K between D 2 and H 2 0, and 8.2 for T 2 
with respect to H 2 32 . For D 2 the melting enthalpy is found 
to be from experiments about 0.07 kcal/mol higher than that 
of water, whereas TIP4PQ/2005 predicts an increase of 0.09 
kcal/mol and TIP4PQ_D20 predicts an increase of about 0.04 
kcal/mol. 

3.3 Isotope effects on the temperature of maxi- 
mum density (TMD) 

Experimentally deuteration of water shifts the TMD by 7.2 
K— , and tritiation by 9.4 K 34 . Our previous results 35 indi- 
cate that deuteration and tritiation of the TIP4PQ/2005 model 
tended to overestimate this shift. Our new models now slightly 
underestimate this shift (see Table [4] and figure In view of 
the fact that the error bar for the TMD is fairly large (±2K), 
these results are reasonable. Especially when one bears in mind 



that an isotopic shift in the TMD is not always present in a num- 
ber of recent model s 32 i 36 , as is the case of the q-TIP4P/F 12 and 
the TTM2.1-F models. 

3.4 The temperature difference between the 
melting point and the TMD 

Of particular interest is the difference between the TMD and 
T m . Experimentally this difference is 3.98 K for H 2 0, 7.5 K 
for D 2 0, and 8.9 K for T 2 0. From our simulations we obtain 
25.7 K for TIP4PQ/2005, 26.8 K for TIP4PQ_D20 and 27.5 
K for TIP4PQ_T20. As can be seen all the models presented 
in this work are unable to describe the difference between the 
temperature of the TMD and the melting point temperature. 
Thus the inclusion of nuclear quantum effects does not solve 
the disagreement with experiment indicating that the origin of 
this failure is the approximate character of the PES. The same is 
true for the q-TIP4P/F which predicts a difference between the 
TMD and the melting point of 26 K (the model has the TMD at 
277 K and the melting point at 25 1 K). For the TTM2. 1 -F mod- 
els^, the difference between the TMD and the melting point is 
even higher (45 K) (the melting point is located at 228 K— 
and the TMD at 273 K^Z) . One can conclude that no model 
designed for path integral simulations thus far is able to repro- 
duce the difference between the TMD and the melting point 
found experimentally. 

3.5 Isotope effects on molar volumes 

Bridgman described that the "molar volume of D2O is always 
greater than that of HjO at the same pressure and tempera- 
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3.6 Isotope effects on the structure 



3 RESULTS 



Model 


AH(T m ) (kcal/mol) 


Pih (T,») 


Pliquid (T m ) 


H2O (experiment) 


1.436 


0.917 


0.999 


D2O (experiment) 


1.509 


1.018 


1.105 


T2O (experiment) 








TIP4PQ/2005 


1.099 


0.919 


0.988 


TIP4PQ/2005 (D 2 0) 


1.189 


1.028 


1.103 


TIP4PQ/2005 (T 2 0) 


1.285 


1.134 


1.212 


TIP4PQ_D20 


1.133 


1.024 


1.091 


TIP4PQ_T20 


1.192 


1.128 


1.214 



Table 5 The change in enthalpy at the melting points of the models along with densities in units of g/cm 3 (experimental values from 
IAPWS-95/NIST Standard Reference Data). 
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Figure 3 Plot of the O-O radial distribution function for water at 290 
K and p = 1 bar. Key: TIP4PQ/2005 red line, TIP4PQ/2005 (D 2 0) 
dashed green line, TIP4PQ/2005 (T 2 0) dashed blue line, 
TIP4PQ _D20 dotted pink line and TIP4PQ.T20 dot-dashed cyan 
line. 



ture"—. Experimentally for ice 1^ this difference was seen to 
be of the order of 0.2% at 220K 40 ' 41 . From our simulations 
of ice I h we obtained 32.410 AVmolecule for TIP4PQ/2005 at 
220K, and 32.303 A 3 /molecule for TIP4PQJD20, also at 220K, 
which is similar to experimental results of 32.367 and 32.429 
for H 2 and D 2 respectively 42 , also at 220K. From this one 
can see that the models used here are unable to capture this (al- 
beit subtle) effect. However, recent ab initio density functional 
theory calculations have been able to reproduce this effect—, 
although there is a w 4% error in the densities themselves. It 
would be interesting to see whether ab initio density functional 
theory calculations are also capable of reproducing the isotopic 
shifts found in the T m and the TMD. 



1.8 
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Figure 4 Plot of the O-H radial distribution function for water at 290 
K and p = 1 bar. The same key as in Fig. [5] 



3.6 Isotope effects on the structure 

In Figures[3]|5]we provide the atom-atom distribution functions 
for O-O, H-H and O-H for each of the models studied. From 
these plots, whose salient features are compiled in Tabled one 
can observe that the structure of the new models, TIP4PQJ320 
and TIP4PQ_T20 is very similar to that of TIP4PQ/2005, more 
so than that of TIP4PQ/2005 (D 2 0) and TIP4PQ/2005 (T 2 0). 
This implies that the structure of the fluid phase of isotopically 
substituted water is almost indistinguishable from that of H 2 0, 
an assumption that experimentalists often make, and one that 
seems to be justified by our simulation results. 

3.7 Isotope effects on the diffusion coefficient 

In the work of Habershon et al.— a H 2 0/D 2 diffusion co- 
efficient ratio of 1.28 was found for the q-TIP4P/F model at 
298K (experimentally it is 1.30). Although it is not possi- 
ble to directly compute the diffusion coefficient from Monte 
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3.8 Heat capacity C ; 



Table 6 O-O, O-H and H-H radial distribution function of water for the various models at 290 K and p = 1 bar. 
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Figure 5 Plot of the H-H radial distribution function for water at 290 
K and p = 1 bar. The same key as in Fig. [5] 



Carlo runs, a rough estimate can be obtained by calculating the 
mean square displacement of the molecules after a fixed num- 
ber of Monte Carlo cycles. In our simulations the ratio of the 
mean square displacement (after 200,000 MC cycles) between 
TIP4PQ/2005 and TIP4PQ/2005 (D 2 0), was 1.33. When we 
compare TIP4PQ/2005 with TIP4PQ_D20 we obtain a ratio of 
1.17, indicating that the new model decreases the differences 
between D 2 and H 2 0, in line with the results for the radial 
distribution function. 



3.8 Heat capacity C p 

The isobaric heat capacity was obtained from a differential 
of the enthalpy with respect to temperature. At 280K the 
values we obtained were (in cal/mol/K) TIP4PQ/2005 17.4, 
TIP4PQ_D20 18.7 and for TIP4PQ.T20 19.5 These results 
compare favourably with the experimental values; H 2 is 
17.9i^andD 2 Ois20.3ii. 

4 Conclusions 

We have seen the "competing quantum effects" interpretation 
of Habershon et al. in action; the TIP4PQ/2005, TIP4PQ/2005 
(D 2 0) and TIP4PQ/2005 (T 2 0) models all have the same ge- 
ometry and charge distribution, thus they all have the same 
dipole moment. When one examines the radial distribu- 
tion functions one can see that the TIP4PQ/2005 (D 2 0) and 
TIP4PQ/2005 (T 2 0) models have stronger features than the 
TIP4PQJ320 and TIP4PQ_T20 models, whose dipole mo- 
ments are smaller. The new models presented here for D 2 
and T 2 were designed by shortening the O-H bond length in 
line with the values presented by Zeidler et al.— . It is worth 
noting that a bond length reduction by as much as 4%, as sug- 
gested by Soper and Benmore^, would probably have led to 
a significant error in our evaluation the isotopic influence on 
the melting point. We have seen that the models studied in this 
work underestimate the melting point by ks 1% and the melt- 
ing enthalpy by rs 30%. This is almost certainly due to the 
approximate nature of the empirical potentials used here, fail- 
ing to reproduce the experimental PES. The situation is more 



7 



REFERENCES 



REFERENCES 



favourable when one considers isotopic shifts. In general we 
have seen that the models qualitatively reproduce the experi- 
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